TOPICS IN MITIGATING RADAR BIAS 
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Abstract 

In this paper, we investigate two topics related to mitigating the effect of radar bias in ballistic missile tracking 
applications. We determine the absolute bias between two radars in polar coordinates when their relative bias is 
given in rectangular coordinates. Using this result, we then obtain the optimized steady-state filter to handle the 
random bias. 

1 Introduction 

There are several facets to the problem of tracking ballistic missiles with radar that require enhanced error correction 
to effectively track threats. In this paper, we obtain the exact form of the bias error for the coordinate transformation 
problem. This result is useful in Ballistic Missile Defense bistatic applications where one sensor is used for launching 
an interceptor, while another is used to track the threat. Thus, the problem of translation between internal sensor 
coordinate frames to a common frame (that is, used by all sensors) is important. The coordinate transformation 
problem from Cartesian to spherical coordinates introduces a bias that, if accounted for, can be corrected in the 
design of a filter. This problem occurs when one has multiple launch platforms, because each local track must be 
formatted for a common reference frame. When bias correction is accomplished correctly, one can improve tracking 
performance of the filter and increase the likelihood that an interceptor can successfully engage a threat. 

2 An Optimized Method of Obtaining Absolute Bias 

Although relative bias calculation can be used to provide correct association of tracks from two sensors, the calculation 
of the absolute bias is required to correct the track state and is needed for track fusion and for producing a Single 
Integrated Air Picture. Methods for obtaining the relative bias between two radars tracking the same ballistic missile 
are presented in Levedahl [2] and Brown, Weisman and Brock [3]. The methods presented in these reports have to 
do with maximizing a likelihood function. The relative biases obtained in these papers are determined in rectangular 
coordinates. In this paper, the absolute bias for the two sensors is calculated from the relative bias by solving a 
minimization problem. The problem is set up to minimize the weighted sum of the two absolute biases while viewing 
the given relative bias as a constraint. 

A point in 3-dimensional space in both rectangular and spherical coordinates 1 is denoted by: 



.r 
V 



and 



r 



, respectively. 



(1) 



The transformations between the coordinates are p = f(n) and tt = f l {p), which are given by: 

r 1 {-y) = 



r cos 9 cos ip 
r cos 6 sin tp 
r sin# 



arctan 



\/ x 2 + y 2 + z 2 
arctan (y/x) 

(z/yjx 2 +y 2 



(2) 



1 Denote yaw (azimuth) by tp, pitch (elevation) by 9. <j> is normally reserved for roll; however, roll is not used here. 
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We need the following definitions: 



Pi 


Target position as seen by sensor f 


P2 


Target position as seen by sensor 2 


Pt 


True target position (unknown) 




Sensor 1 bias 


B 2 


Sensor 2 bias 


Br 


Relative bias 


P1T02 


Sensor 2 position from sensor 1 



Pl.ENU(l) — { x l, 2/1, z i)'eNU(1)' Pl,ENU(2) — {x 2 , 2/2, Z 2)' ENU ( 2 ) ■ Bl,ENU(l) ~ {^- x l , ^2/1 , ^ z l)' E NU(l) ' B 2,ENU(2) = 

(Ax 2 , Aj/2, Az2)' ENU ( 2 )- (ENU denotes the East North Up coordinate system.) Thus, we have in the sensor coordi- 



nates 

Pt,ENU(1) — Pl,ENU(l) + Bl,ENU{l) 
Pt,ENU(2) = P2,ENU(2) + B2,ENU(2) ■ 

If we use an ENU coordinate system located at sensor 1, (4) becomes 

Pt,ENU(1) = PlT02,ENU(l) + ^2,_E7V(7(1) + B 2 ^ENU(1) 



(3) 
(4) 

(5) 



where P\to2.enu(\) is the position vector from the first sensor to the second sensor in ENU{1). The relative bias 
in ENU(l) is 

Br,enu{i) — B 2 ^enu(i) ~ B 1ENU ( 1 ) 

P2,ENU(l)) ~ {Pt,ENU(1) ~ Pl,ENU(l)) 



(Pt,ENU(1) ~ PlT02,ENU(l) 



Pi. 



EN (7(1) 



P, 



1T02,ENU(1) 



2,ENU(1) 



(6) 



We consider the coordinate transformations to allow us to go from ENU to radar-face coordinates for a particular 
sensor. Each sensor has its own face and ENU coordinate systems. The face coordinate system (denoted FACE) of 
a sensor is related to the ENU coordinate system of a sensor by the following transformation: 



ENU(i)2FACE(i) — 



cos 9i 


— sin 9; 



cos 9i cos f/'i 

— sin-0^ 
— sin 9i cosV'i 



sin 9i 

1 

cos6»i 



cos 9i sin ip i 

COST/'i 

— sin 9i sin ip i 



cos^j sinf/'i 
— sin if) i cos ipi 


sin 8i 


cos 9, 



(7) 



where i = 1,2. We also have that 



l FACE(i)2ENU(i) 



cos 9i cos ipi — sin tp i — sin 9i cos 



cos 6i sin^j 
sin 9j 



cos ^ 




sin 9i sin ip i 
cos 9j 



(8) 



which is the transpose of (7). We can also have the matrix T ENU ^ 2 FACE(j), which is 

cos 6i j cos 4 cos 9i.i sin ?/>,• sm9 



l ENU(i)2FACE(j) 



- sin i/j^j 
sin^i.j cos ip i : j 



cos 9ij sin 
cos tp^- 
— sin 9ij sin ^ j 







The absolute (as opposed to relative) bias can be expressed in the face coordinates: 

Bi^FACE(i) = Ar • T? r + Ac A ■ u^X + Ac B ■ 



(9) 



(10) 
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where is the unit vector in the range coordinate and u c X, u c b are the two cross range coordinate unit vectors. 
Substituting prAip = Aca and ptA8 = Acb where pxi = ||-Pr(*)ll ( see note 2 ), the distance from sensor to the 
target, we get 

Buface(i) = Ar, L ■ u^+pTiAipt ■ u^X+p Tl A9 l ■ (11) 



D 



i,ENU(i) — 



COS Oi COS 

cos 9i sin 



— sin — sin Oi cos rp i 
cos — sin Oi sin tp i 
cos 6; 



An 

PTtAlpt 

PnAOi 



Ari cos Oi cos ip i — Atp i (sin ip i ) • pn — AOi (sin Oi cos ip i ) ■ pTi 
Ari cos Oi sin 'ip i + (cos ^ i ) • pri — AOi (sin #j sin ^ ) • p Ti 
An sin Oi + AOi (cos 6>; ) • p T i 



(12) 



We can obviously obtain B iENU ^ (for i not necessarily equal to j) if needed. The quantities Ari, A^ 1 , A#i, Ar 2 , 
A^ 2) A6>2 are the ones we minimize. To refer to these as a group we on occasion write e = (An, Atp 1 , AOi, Ar 2 , Aip 2 , A9 2 ). 
We need tolerances or costs for the sensor biases. These are expressed in spherical coordinates. 



k r i 


Sensor 1 range bias cost, unitless 


k r2 


Sensor 2 range bias cost, unitless 




Sensor 1 azimuth bias cost, meters 


k^2 


Sensor 2 azimuth bias cost, meters 


kei 


Sensor 1 elevation bias cost, meters 


ke2 


Sensor 2 elevation bias cost, meters 



2.1 Problem Statement 

We want to compute the minimum (absolute) bias cost for the two sensors when there are known (computed) 
expressions for the relative bias. The given relative bias is expressed in ENU rectangular coordinates. We compute 
the minimum absolute bias in spherical coordinates. The relative bias in rectangular coordinates contrasted with the 
absolute bias in spherical coordinates allows us to formulate this as a minimization problem. We view the relative 
bias as a constraint. We use a quadratic cost: 



F = f ■ (An) 2 + \ • (A^f + \ • (AO.f + f • (Ar 2 ) 2 + \ ■ (A0 2 ) 2 + f • [A9 2 f 



(13) 



So that the addition in (13) is permissible, we have that k ri , k r2 are unitless and k$ 1 , k$ 1 , k^, 2 , kg 2 are in meters. 
We note F may be rewritten in the form 





r 2/k* n 








-l 


An 


F = [ An Aip 1 A6»i ] 












Atpx 








Vkl _ 




A0 1 



+ [ Ar 2 A-0 2 A6 2 



Vk 2 r2 






o 






Vkl 



-, -1 


Ar 2 




A^ 2 




A0 2 



(14) 



which we recognize as being in the form of a Mahalanobis distance. We note that the Mahalanobis distance comes 
up in the Levedahl/Lincoln Labs work ([2] and [3]) when the log is taken of the Gaussian distribution. The cost F 
is minimized subject to this equality constraint: 



G{B) 



IB. 



2,ENU(1) 



B 



1,ENU(1) j 



B 



R,ENU(1) 



= 



(15) 



Thus, we have 



G(B) 



Ar 2 cos #2 cosf/' 2 — Aip 2 (sin^ 2 ) • PT2 — A0 2 (sin #2 cosf/' 2 ) • PT2 
Ar 2 cos O2 sin tp 2 + Aip 2 (cos ip 2 ) • pr2 — A9 2 (sin 2 sin tp 2 ) • pr2 
Ar 2 sin 2 + A9 2 (cos 9 2 ) • p T 2 



2 True position is not available. When applying this method measured position is used for this calculation instead. 
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Ari cos 6*i cosV'i — A0i (sini/^) • p^i — A0i (sin0i cos ^i) ■ Pti 
Ari cos 0i sinV'i + A0j (cos0>i) • Pti — A0i (sin0i sin0i) • pn 
An sin^i + AO i (cos0i) -pri 



-Br = 



(16) 



where all of the terms in (16) reside entirely in one or the other of the two ENU coordinate systems. We see that 
G(B) gives that the difference between the two absolute biases (whatever they may be) is equal to the relative bias. 
Also, we note that (16) is affinc. Another equivalent representation for G(B) is 



G{B) = A(p T2 ,0 2 ,0 2 ) 



Ar 2 

A0 2 
A0 2 



A(pri,^i,0l) 



Ari 
AVi 
A6»i 



Br 



where 



A(pri,Vi,(9i) 



COS 0\ COS'i/'i 

cos 0\ sin^ 
sin^i 



-sin^i - Pti 
cosip 1 • pti 




cos 2 cos tp 2 
cos0 2 sin -02 
sin 2 



- sin 2 ■ p T2 

COS 2 • p T2 





A(p T2 ,ll> 2 ,0 2 ) - 

Setting G(B) = 0, we solve for Ar 2 , A0 2 , A0 2 

■ A' 1 ( PT2 ,ip 2 ,0 2 ) I A {pti, 0i, 0i) 



- sin 0\ cos 0i • pti 

- sin#i sin^ • pxi 

cos0i ■ pti 

sin #2 cos 02 • Pti 
sin 6*2 sin 2 • Pti 
cos 2 • Pti 



(17) 



(18) 



(19) 



Ar 2 

A0 2 

A0 2 



Ari 
A0i 
A0i 



(20) 



(provided j»t2 0.) This vector equality constraint (16) can be written in the form of three scalar equality constraints 
Ge(B) = Ar 2 cos 2 cos 2 — A0 2 sin 2 ' PT2 — A0 2 sin 2 cos 2 1 PT2 

—Ari cos 0i cos 0-l + Aip 1 sin X ■ pti + A6\ sin 6\ cos 0^ • P ti — Br E (21) 
Gat (B) — Ar 2 cos 9 2 sin 2 + A0 2 cos 2 ' PT2 — A9 2 sin 9 2 sin 2 ' PT2 

—Ari cos 0i sin^ — Atp 1 cos0i • pri + A0i sin0i sin0 1 • P ti — Brat (22) 
Gu{B) = Ar 2 sin6> 2 + A0 2 cos0 2 -PT2 - Ari sin0i - A9 1 cos0 1 ■ p T1 - B RU . (23) 

2.2 Solving The Minimization Problem 

To solve this minimization problem, we need to take a few derivatives. We need the gradient of the function to be 
minimized. We also need the gradient of the constraint, which is an equality constraint in this case. 



dF/dAn 




' k 2 


•Ari 


dF/dAip 1 




k 2 


AVi 


dF/dA9 1 






A0i 


dF/dAr 2 




k 2 
K r 2 


• Ar 2 


dF/dAip 2 
dF/dA9 2 




k 2 


A0 2 




- ^02 


A0 2 



(24) 



dG E /dAn 
dG E /dA^p 1 
dG E /dA9 1 
dG E /dAr 2 
dG E /dAiP 2 
dG E /dA0 2 



— COS 01 COS 0^ 

sin t/0 -p T1 
sin 0i cos 0i • pti 
cos 02 cos 02 
- sin 02 • PT2 
— sin 02 cos 2 ' PT2 



(25) 
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VG 



N 



dG N /dAn 
dGN/dA^ 
dG N /dA9 1 
dG N /dAr 2 
8Gn / dA%p 2 
dG N /dA9 2 

dGu/dAn 
dG u /dAip 1 
dG u /dA9 1 
dGu/dAr 2 
dGu/dAxP 2 
dG u /dA0 2 

We are looking for an optimal solution located at the point e* 
Kuhn- Tucker conditions that stipulate the optimal solution e 
there exist numbers a\,a 2 , Og such that 



(26) 



VGu = 



— cos 6*i sin^ 
-cos?/>i -Pti 
sin 0i sin0i • pn 
cos 9 2 sin -02 
cos 2 ■ p T2 
- sin 9 2 sin V> 2 • pt 2 

— sin 


— pn • cos 9\ 
sin 2 


p T2 ■ COS 2 

= ( Ar I , A0* , A0* , Ar* , A^* 2 , A9* 2 ) . We employ the 
should satisfy these equality constraints for e and 



(27) 



VF (e*) = a\ ■ VG E (e*) + • VGjy (e*) + a* 3 ■ VG V (e*) 



(28) 



The gradients VG^ (e*) , VGat (e*) , VG[/ (e*) are linearly independent. Taking an inventory of the equations and 
unknowns, we see that there are 9 unknowns (e, 01,02,03) and 9 equations (3 from the equality constraint and 6 
from the above equation). We may be able to find the solution. Since the cost F is quadratic and the constraint G 
is afhne, the necessary conditions we give for optimality are also sufficient conditions and an optimal solution e* is 
a global optimal solution. Equation (28) in longhand is: 



r k 2 


An 




— COS 01 COSf/'i 




— cos 0i sin^i 


k 2 

>i 






sin 0i - pn 




-cos^i - Pti 


k 2 


A0i 




sin 0i cost/ 1 ! • Pti 


+ a 2 


sin 0i simp 1 ■ pri 


k 2 

K r 2 


Ar 2 


= ai 


COS 02 COS 2 


cos 02 sin 2 


k % 2 


A0 2 




- sin ip 2 ■ p T2 




COS 2 ■ PT2 


- ^02 


A0 2 




— sin 2 cos ip 2 ■ pr 2 




— sin 2 sin ip 2 ■ pt 2 



03 



— sin 0i 


— pti ■ cos 0i 
sin 2 


PT2 ■ COS 02 

The right hand side of (29) may be written in the form of the product of two matrices Mi and M 2 , 



(29) 



— cos 01 cos tp 1 
sini/>i -p T i 

sin 0i cos -01 • pti 
cos 02 cos 2 

- sin V> 2 • PT2 

— sin 02 cos ip 2 ■ pti 



— cos 0i sin0i 
-cos^i - pti 
sin 0i sin ^i ■ Pti 
cos 02 sin 2 
cos tp 2 ■ pt 2 
— sin 02 sin -02 ■ PT2 



— sin 0i 


— pti ■ cos 0i 
sin 2 


PT2 ■ COS 2 







"Mi" 




Oi 


a 2 




M 2 _ 




2 






. fl 3 . 








. fl 3 . 



(30) 



where 



— cos 0i cos 0i 
M\ = sin -01 • pti 

sin 0i cos 0i • pti 

cos 02 cos 2 
M 2 = - sin0 2 ' PT2 

- sin 02 cos 02 ' PT2 

Note that oi, a 2 , 03 have the units of meters. Let 

Di = 



— cos 0i sin 0i 
-cos 0i - pti 
sin 0i sintp 1 ■ pti 

cos 02 sin 02 

COS 02 • PT2 

— sin 2 sin 2 ' PT2 PT2 • cos 2 



— sin 0i 


~Pti ■ cos 01 

sin 2 




D 2 = 



' k 2 



k 2 





k 2 












(31) 
(32) 

(33) 
(34) 
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Rewriting the left hand side of (29), 



" k 2 
















k 2 

















k 2 

K 6i 

















k 2 

K r 2 

















k 2 






















Mi 
M 2 



An 

A6»i 
Ar 2 

AV> 2 
A0 2 

ai 

02 
03 



3 .3 ^2 



An 

A6»i 

Ar 2 

AV 2 
A6 2 



Hence, we have 





An 




ai 






= Mi 


02 




. A6, i . 








Ar 2 




ai 


D 2 


AV-2 


= M 2 


02 




A9 2 




03 



or, 



Ar 2 
AV> 2 
A0 2 



= D 2 X M 2 M^D X 



An 
AV'i 
A6»i 



Substituting (38) into (20) yields 
D 2 1 M 2 M^ 1 D 1 



An 
A6>! 



= ^ 1 (p T 2,il> 2 ,0 2 ) A(pri,Vi^i) 



An 

A^i 
A6»i 



+ 5/ 



so we get 



(D^ 1 M 2 M^ 1 D 1 - A- 1 {vti^M A (pruipM) 



An 
AVi 
A6»i 



A 1 (p T 2,ip 2 ,0 2 )B R 



An 
AV'i 
A6»i 



= {D^ 1 M 2 M^ 1 D 1 -A- 1 ( P T2,il>MA( P Ti,il> 1 ,0i)) 1 ' A- l {pT2,tl>MB R , 
which allows us to obtain (An, AV>i, A#i). Finally, substituting (41) into (38) we get (Ar 2 , Atp 2 , A9 2 ). 



2.3 Numerical Examples 

We illustrate this idea with a few examples. 

2.3.1 a. 



INPUT 


OUTPUT 


B R = [ 200 500 300]' 
Pti = 25000 

Vh = o 

6»i = 7.8540e-001 
PT2 = 50000 

V 2 = 

9 2 = 2.3562c+000 


Cost = 1.6250e+004 
An = -1.7678e+002 
A^i = -1.0000e-002 
A6»i = -1.4142e-003 
Ar 2 = 3.5355c+001 
AV-2 = 5.0000e-003 
A6>2 = -3.5355e-003 



INPUT CONTINUED 


OUTPUT 


fcl = 1.2500c+009 = 2 * PT1 2 . 
fcfi = 1.2500e+009 

fc2 = 5.0000c+009 = 2 * PT2 2 . 
kj 2 = 5.0000c+009 





2.3.2 b. 



INPUT Same as with a. but with 


OUTPUT 


B R = [ 200 500]' 


Cost = 3.6250e+004 
An = -2.4749c+002 

AV>! = 

A6t = -4.2426e-003 
Ar 2 = 1.0607e+002 

Aip 2 = 

A6> 2 = -4.9497e-003 



2.3.3 c. 



INPUT Same as with a. but with 


OUTPUT 


tp 2 = 7T 

2 = tt/4 


Cost = 1.6250c+004 
Ari = -1.7678c+002 
Aijj 1 = -l.OOOOe-002 
A6»i = -1.4142c-003 
Ar 2 = 3.5355e+001 
Ai/> 2 = -5.0000e-003 
A6>2 = 3.5355e-003 



The input for this case is a variation of the input in a. Note that the output is the same as with a except for a 
sign swap between Aip 2 and A9 2 to account for the orientation difference of the "2" coordinates. 
2.3.4 d. 



INPUT Same as with a. but with 


OUTPUT 


2 = tt/4 


Cost = 5.5625e+004 
Ari = -1.7678c+002 
AV>! = -l.OOOOc-002 
A0i = -1.4142c-003 
Ar 2 = 2.8284c+002 
Aip 2 = -2.0000c-003 
A6> 2 = -1.4142c-003 



3 An Optimized Reduced-State Filter For Unknown Bias 

A novel technique for calculating a steady-state reduced-order filter to track a maneuvering target is presented by 
Mookerjee and Reifler [4]. The filter they derive is optimized for performance with a stochastic acceleration. In this 
paper, this technique is modified to derive a steady-state filter that is optimized for performance with a stochastic 
measurement bias. Similar to [4], the filter developed here is a reduced-state filter. We can see what a reduced-state 
filter is by considering [6] and [5]. In these reports, we estimate the position and velocity of an aircraft (a Beechcraft 
1900) with DMEs (distance measuring equipment), an INS (inertial navigation system) and a barometric altimeter. 
The filter (in [6]) and the smoother (in [5]) were designed with a state-to-estimate range bias in each DME (up to 5 
were used), a state-to-estimate INS drift, and a state-to-estimate bias in the baro. The filter (or smoother) ran with 
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these additional bias states in tow (i.e., in addition to the position and velocity states). (The results in [6] and [5] 
achieved the design goals in position and velocity accuracy.) 

It appears likely the design goals of [4] and this paper are competing design goals. The design methods discussed 
in [7], which are based on the Bode gain-phase relationship, possibly could be brought to bear to quantify a possible 
trade-off on the design goals of this paper and [4]. We don't cover such trade-offs in this paper, but it could be a 
problem for future investigations. The classical control concepts of the sensitivity function and the complimentary 
sensitivity function come to mind. 

We use discrete time dynamical equations. It is fair to consider our state and output (dynamical) equations to 
be the dual (in the control theory sense) of the state and output equations, Equations (8) 3 and (5). Compared to 
the dynamical equations in [4] , we eliminate the unknown acceleration from the state equation and add an unknown 
bias in the output (measurement) equation, the typical dual situation. We have: 

x(k + 1) = <f>x(k) +m(k) (42) 

z (jfe) = H ■ x (jfe) + n (jfe) + Wu (x (jfe) , A) . (43) 

The state x(k) at time k is of dimension n and the state transition matrix $ is of dimension n by n. The output 
z{k) at time k is of dimension q and the output matrix H is of dimension q by n. The process noise term m(k) is of 
dimension n with covariance Q. The measurement noise term n(k) is of dimension q with covariance N. The bias 
matrix W is q by to. The bias function u is JJ™ x W —* 5i m , and we have that the bias A is a p-dimensional random 
vector with mean A and covariance A. 

The time update equation, using (42), is simply 

x(k + l\k) = <5>x(k\k) . (44) 

The measurement update equation becomes 

x(k + l\k + l) = x(k + l\k) + K (z(k + l) - 

where K is the n by q measurement, or Kalman, gain matrix 
position gain a and velocity gain (3 substitute for K. 

3.1 Filter Development - General Case 

In this subsection, we develop the filter equations for the general case. The development in this section is (basically) 
dual (dual in the sense of control theory) to Section III in [4] . The error is defined as (we develop the errors analogous 
to (27) and (32)): 

e(fc + l|fc+l) =x{k + l)-x{k + l\k + l) (46) 
= x (k + 1) - x (jfe + l|jfe) - K (z (jfe + 1) - Hx(k + l\k)- Wu (jc (jfe + l|jfe) , A)) 

= x(k + l)-x(k + l\k) 
-K (Hx (jfe + 1) + n (jfe + 1) + Wu (x (jfe + 1) , A) - Hx (jfe + l|jfe) - Wu (x(k + l|jfe) , A)) . 

Continuing, 

e {k + l\k + 1) = x (k + 1) - KHx (k + 1) - Kn (k + 1)- KWu (x (k + 1), A) 
-x (k + l\k) + KHx(k + l|jfe) + KWu (x (jfe + l|fc) , A) 
= <I>x (k) + to (k) - KB (<f>x (k) + m (fc)) - Kn (k + 1) - KWu ($x (k) , A) 
-$x(k\k) + KH$x(k\k) + KWu ($x(k\k) , A) 
= (I - KB) $ [x (k) - x (k\k)) + KB) to (jfe) -Kn(k + 1) 
-KW(u($x{k),\)-u($x(k\k),\)) . 

So we have 

e (jfe + l\k + 1) = L$e (k\k) + Lm (k) - K (l^Au fe | fc + n (k + 1)) ; (47) 

3 In this section, we often refer to equations from [4]. Hence we adopt the convention that all equation references appearing in bold 
typeface are to equations in [4] . 



- Hx (k + l\k) - Wu (x (k + l\k) , (45) 
. In the steady-state case, which is discussed below, the 
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where 

an n by n matrix, and 



L=(I- KH) 



Au k \ k = u (<&x (k) , A) - u ($x (k\k) , A) . 
We can make the linear approximation 



du 

Aufe i fc " Tx 



$Ax + - 



AA 



so 
and 

We obtain the result: 



x— x{k\k) .X— X 

Ax = e (k\k) = x (k) - x (k\k) 
AA = A - A . 



x—x(k\k),X— X 



(48) 
(49) 

(50) 

(51) 
(52) 



- KW \Tx 



e (k + l\k + 1) = L<S>e (fc|fc) + Lm (k) 

AA) -Kn(k + 1) 



*s(k\k)+^ 



L-KW 



du 
dx 



x— x(k\k),X= A 

du 



x—x(k\k),X—X 



dX 



where 



an n by n matrix, and 



x—x(k\k)^X— A 

^ $e(fc|jfe) +Lm (k)-KW 
= Fe (k\k) + Lm (k) + CAA - Kn (k + 1) 



AA - Kn (fc+1) 



F = \ L — KW 



C = -KW 



dx 

du 
dX 



x— x(k\k),X— X . 



x— x(k\k),X— X 



(53) 
(54) 

(55) 



an n by p matrix. 

We now implement the observation made in [4] that the error e (k\k) may be viewed as consisting of two compo- 
nents. The first component of error, , is due to the process noise m and the measurement noise n. The second 
component of error, e^ 2 \ is due to the measurement bias. To the extent that the linear approximation is valid, a linear 
analysis holds. That is, the two error inputs may be treated in separate equations by applying the superposition 
principle of linear analysis. 



e (1) (k + l\k + 1) = (k\k) + Lm (k) -Kn(k + 1) 

(k + l\k + 1) = FeM (k\k) + C • AA . 



(56) 
(57) 



These equations are comparable to (33) and (34). 

In addition, we require update equations for the total covariance and the covariance of (k\k). Using (42), the 
time update equation for (k\k) is 



M (k + l\k) = E s w (k + l|fc) £ W (k + l\k)' 
= $M(fc|/c)$' + Q . 

From (56) the combined (measurement and time) update for the covariance of (k\k) is 



(58) 



M (k + l\k + 1) = E \e (1) (k + l\k + 1) gW (k + l\k + 1)' 



9 



E 



(Ve (1) (k\k) + Lm (k) -Kn(k + lfj (k\k) + Lm (k) -Kn(k + 1)) ' 

and we use E [n (k) n (/)'] = for k ^ I giving E [Fe« {k\k) (Kn (k + 1))'] = 0. Hence, 

M (k + l\k + 1) = E \Fe {1) (k\k) e w {k\k)' F' + Lm (k) m (k)' L' + Kn(k+l)n(k + 1)' K ' 

= FM (k\k) F' + LQL' + KNK' . (59) 

Working towards update equations for the total covariance, we firstly define the n by p matrices D (k\k) and 
D(k + l\k) as 

e {2) (k\k) = D (k\k) ■ AA . (60) 

We can define D (k\k) in this way since in view of our linearized analysis, the system output (e' 2 ) (k\k)) is a linear 
function of the system input (AA). We proceed by defining 



D(k + l\k) = FD(k\k) . 



(61) 



In (60), and AA are known quantities (the equation defines D (k\k)). In (61), F and D (k\k) are known quantities. 
Then, substituting (60) into (57), we obtain 



D (k + l\k + 1) • AA = FD (k\k) ■ AA + C ■ AA 

= D(k + l|fc) • AA + C- AA , 
and subsequently (assuming (62) holds for all AA.) 

D(k + l\k+l) = D(k + l\k) + C . 



(62) 



(63) 



Let S be the total error covariance. By superposition, we get the total error by the addition of the two error 
terms. We also make the observation that since the two errors, and e^ 2 \ originate from independent sources 
they remain independent for all times k. Looking at S*, 



= E 



S (k + l\k) = E [e (k + l\k) e (k + l\k)'] 
(k + l\k) + e {2) {k + l|fc)) (e (1) (k + l\k) + e {2) (k + l|fc))' 



= E e«(fc+ l|fc)e«(fc+l|fc)' +E e (2) (k + l\k) e (2) (k + l\k)' 

= M(k + l\k) + E \e {2) (k + l\k) e (2) (k + l\k)' 

= M (k + l\k) + E [$D (k\k) AA • AX'D (k\k)' , 
using (42), (44) and (60). Hence, 

S (k + l\k) = M (k + l\k) + $L> (k\k) E [AAAA'] D (k\k)' . 



Finally, 



S{k + l\k) = M (k + l\k) + $D (k\k) AD (k\k)' 



(64) 



Basically, this is the same result as (19). 

We next obtain the measurement update for S: 

S(k + l\k + l) = E [e(k+ l\k+ l)s(k+ l\k + 1)'] 

= E[(x(k + l)-x(k + l\k + 1)) {x {k + 1) - x (k + l\k + 1))'] 
= E[(x(k + l)-x(k + l\k) -K[z(k + 1)-Hx(k+ l\k) - Wu (x(k + l\k) , A)]) (ditto)'] 
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= E[(x(k + l)-x(k + l\k) 
-K [Hx (k + 1) + n (k + 1) + Wu (x (k + 1) , A) - Hx (k + l\k) - Wu (x (k + l\k) , A)]) (ditto)'] 
= E [{x (k + 1) - x(k + l\k) - KH [x (k + 1) - x {k + l\k)} 
-K [n {k + 1) + Wu (x (k + 1) , A) - Wu (x (k + l\k) , A)]) (ditto)'] 

= E [{x (k + 1) - x (k + l\k) - KH [x (k + 1) - x (k + l\k)} 
-K [n (k + 1) + W (u {x (k + 1) , A) - u (x (k + l\k) , A))] } {ditto}'] 
= E [{x (k + 1) - x(k + l\k) - KH (x(k + l)-x(k + l\k)) 
-K [n (k + 1) + W (u {$x (k) ,X)-u (<f>x (k\k) , A))] } {ditto}'] 
= E [{x (k + 1) - x (k + l\k) - KH (x {k + 1) - x(k + l\k)) - K [n (k + 1) + W (Au k \ k )] } {ditto}'] 
We take this next step only to the extent of the approximation, 



S(k + l\k + l) = E 



(I - KH) (x(k + l)-x(k+ l\k)) - K 
x {ditto}'] 



„(/, + ]., + ir (g$ e (fc|fc) + ^AA 



E 

= E 



(I - KH) e(k + l\k) - K 
du 



I - KH - KW 



dx 



Let 



and 



Then 



n(k+ 1) + W 
e(k + l\k)-K 

H = H + W 



du . , , . , . du 
-e(k + m + -AX 



rlii 

n(k + l) + W—A\ 



{ditto}' 
{ditto}' 



du 
dx 

n = n + w£a£'w. 

OA OA 



Now 



S (k + l\k + 1) = (i - KH} S(k + l\k) (i - KH^j + KNK' 
- (/ - KH) E[e(k + l\k) AA'] (w^j K' 
-K (w^j E [AXe (k + l\k)'] (/ - KH)' . 

E[e{k + l\k) AX'] = E [(e (1) (k + l\k) + e {2) (k + l|fc)) AA' 



E 



e {2) (fc + l|fc)AA' 



= <PE 



e {2) (k\k) AA' 



Hence, 



= <&E [D (k\k) AAAA'] = <PD (k\k) A . 

S (k + l\k + 1) = (j - KH) S(k + l\k) (i - KH)' + KNK' 
-(l-KH)$D(k\k)A(w^ K' - K (w^j AD (k\k)' ' (i - KH)' . 



(65) 
(66) 



(67) 



Equation (67) is similar in form as (37). We select K so as to minimize the trace of S (k + l\k + 1): tr (S (k + l\k + 1)). 
We minimize tr (S (k + l\k + 1)) because for positive definite matrices trace going to zero implies the matrix L 2 norm 
goes to zero. Let P be a positive definite n by n matrix, we have 



\P\\ < tr(P) < n ||P|| 
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So minimizing tr (P), gives a smaller upper bound for ||P||. We find the optimal K (K that minimizes tr (S (k + l\k + 1))) 
by taking derivatives with respect to K, setting the result to zero and solving for K. We recall the following facts: 
Let A be a matrix independent of K, 

^-tr{KAK') = K{A + A') ; 



dK 



^r(KA)=A'; 
-Jj-j^tr (K'A) = A ; 
^r(AK)=A' ] 
^tr(AK')=A. 



Using these facts on the terms of (67), 



-^tr {KHS{k + l\k)H'K'^ = 2KHS (k + l\k) H' ; 



d 

OK 
d 
~dK 



tr (kNK 1 ^ = 2KN ; 
tr(S(k+l\k)) = ; 



-^tr [S (k + l\k) H'K'} =S(k + l\k) H' 



dK 
d 



—tr (KHS (k + l|fc)) = S (k + l\k) H' ; 



dK 

^Ltr ($D (k\k) A (w^j K^j = <S>D (k\k) A (w^ 
-^-tr^KH^D{k\k)k{w^j K^j =K^H^D{k\k)k{w^j + (w^J AD (k\k)' VH'^J ; 

d . f„f„ T du\ k ^,^^ l ^ l ~ l \_^((Bu\.^ n , 1 . l ^~ l . 



— tr [k (w-q^j ad ( fc l fc ) &h'k'J = K y y W dx) AD ( k \ k > ®' H ' + H ® D A y^ax 

We differentiate the trace of S (k + l\k + 1) as it is represented in (67) by K and set the result equal to zero, 
—tr (S(k + l\k+ 1)) = 2KHS (k + l\k) H' + 2KN + - S (k + l\k) H' - S (k + l\k) H' 

-^D (k\k) A (w^j + K ^H$D (k\k) A (w^j + (w^j AD (k\k)' &H^J 
-$D (k\k) A (w^) +k( (w^) AD (k\k)' + H<S>D (k\k) A (w^) J = . 



After some algebra, 



K \ US(k - i i. ) W + A - U<I>D (k\l, ) A [ + (w^- ) AD (k\k)' 
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S(k- I k)IJ' + ^D(k\k) \ 



(68) 



The optimal filter gain is 



K - \lJS(k- // \ //•'•/' \ i + (w^J ) AD(k\k)'&H' 



> f .S'i7,-+ W - <I>/J (/,-|/,-) A ( 



(69) 



which is an n by g matrix. 



3.2 Filter Development - Steady-State Case 

We next examine the steady-state case of our problem. Referencing equations (41-43) we have 



$ = 



1 T 
1 



H=[l 
W = 1 . 



(70) 

(71) 
(72) 



Hence, with p as position, v as velocity, and z as the measurement, the state transition and output equations for the 
steady-state case are 



p(k + l) 
v(k+l) 



1 T 
1 



p(k) 
v(k) 



+ 



p{k) 
v{k) 



z(k)=[ 1 ] 

Setting u (x, A) = A, our linear approximations from (50) become 

£ -=[° °1 

ux x=x{k\k),\=\ 

and 



m (k) 

+ n(k) + u(x(k),X(k)) . 



du 



i . 



x— x(k\k) ,X— A 



Then, substituting (75) into (65) 

H = H + W d- X = ^ 1 + °] = [ 1 °]= ff 

and substituting (76) into (66) 



N = N + VK— A— W' = JV+l-l-A-l-l = iV + A. 
oA oA 



We have that the steady-state filter gain is 



K = 



a 
(3/T 



(73) 
(74) 

(75) 
(76) 

(77) 

(78) 
(79) 



As mentioned previously, (79) is where a and [3 fit in for the Kalman gain matrix K as given by (69). These gains 
are obtained by computing the steady-state values for all the variables in (69). The major objective of this section 
is to find a relationship between a and (3. 
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1 
1 



The steady-state version of L from (48), L, is 
L=(I- KH) 
The steady-state version of F from (54), F, is 
F = 



a 
/3/T 



[1 0] 



l-a 

-(3/T 1 





'l-a " 




a 


( 


-(3/T 1 




a /t 





'IT' 




"l-a (l-a)T 


) 


1 




-/3/T 1 - /3 



The eigenvalues of F arc 



1 • [ 



Ai, 2 = 1 - ± i^2a/3-4/3 + a 2 +/3 2 



Then, referring to (55), the steady-state version of C, C, is 



C=- 



(3/T 



• 1 • 1 = - 



a 
(3/T 



The measurement updated steady-state covariance M, referring to (59), is 



M = lim M(k\k) = FMF + LQL +KNK 



Superimpose the two noise terms by letting 
and solve: 



M = M Q + M N 

J7\ =fm n f' + knk' 
Mq = fm q f' + LQL' . 



(80) 

(81) 
(82) 

(83) 

(84) 

(85) 

(86) 
(87) 



Comparing (79), (81) and (86) to (46), (47) and (48), we see that our solution for Mn is of the same form as (49). 
Consequently, 



Mat 



N 



2a 2 + 2(3- 3a/3 (3 {2a - (3) /T 
(3 (2a - (3) /T 2(3 2 /T 2 



a (4 - 2a - (3) 

The solution of Mq remains to be determined. Substituting (80), (81) and (73) into (87) gives 



(88) 



M Q = 
In longhand, 



1 — a (1 — a)T 

-a /t 1 - (3 



M Q 



1 - a -(3/T 
(l-a)T 1-/3 





'l-a " 




' 




+ 


-/3/T 1 _ 




q 2 2 





1 - a -(3/T 
1 



Mq 



muQ m 12 Q 
mi2Q m 2 2Q 



(l-a) (muQ + 2Tmi2Q + T 2 m 2 2Q) 
(1 - a) (-(3m 11Q /T + (1 - 2(3) m 12Q + T (1 - (3) m 22Q ) 

(1 - a) {-(3m 11Q lT + (1 - 2(3) rn 12Q + T (1 - (3) m 22Q ) 
(/3 2 /T 2 ) m llQ + {2(3 ((3-1) /T) m 12Q + (l - 2/3 + (3 2 ) m 22Q 





q 22 



(89) 



Hence. 



UiQ - (1 - a) (ranQ + 2Tmi2Q + T 2 to 22 q) 



mi2Q - (1 - a) {-(3m 11Q /T + (1 - 2/3) m 12Q + T (1 - (3) m 22Q ) 

mi 2Q - (1 - a) {-(3m 11Q /T + (1 - 2/3) 7n 12Q + T (1 - (3) m 22Q ) 
"T 22Q - (/3 2 /T 2 ) mug - (2/3 (/3 - 1) /T) m 12Q - (l - 2(3 + (3 2 ) m 22Q 



q 22 



(90) 
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We get three equations in the three unknowns mug, m\ 2 Q and m 22 Q, 



i-(i- a y 



-2(l-a) 2 T -(l-a) 2 T 2 
(l-a)0/T 1 - (1 - a) (1 - 20) - (1 - a) (1 - 0) T 
-{0 2 /T 2 ) -20 {0-1) /T l-(l-20 + 2 ) 




922 



muQ 
m 12 Q 
m 2 2Q 



Taking the matrix inverse to solve for Mq, 



muQ 




' 1 - 






(1 


_ ra 22Q 







-2 (1 - a) 2 T 



2^2 



-1 












. <?22 



The determinant of this matrix, 4a0 — a0 2 — 2a 2 = a0 (4 — — 2a), should not be zero for the inverse to exist. 
This is satisfied by these conditions: 

1. a ^ 

2. 0^0 . (91) 

3. 0^4- 2a 



If the determinant is not zero, we can obtain the solution: 









= 922 ' 


_ m 2 2Q 





In matrix form 



<722 



(-4a0 + a0 2 + 2a 2 0) 
And finally M is obtained from (85), 



T 2 {-2 + ha -4a 2 + a 3 ) 
T (-2a + 0-a0 + 3a 2 - a 3 ) 
(-20 + 2a0 - 2a 2 + a 3 ) 



T 2 (-2 + 5a - 4a 2 + a 3 ) 
T(-2a + 0-a0 + 3a 2 -a 3 ) 



I (-4a0 + a0 2 + 2a 2 0) 



T(-2a + 0-a0 + 3a 2 - a 
(-20 + 2a0 - 2a 2 + a 3 ) 



(92) 



(93) 



and (93): 



M 



N 



(4-2a-0) 



2a 2 + 20- 3a0 (2a - 0) /T 
0(2a-0)/T 20 2 /T 2 



922 



(-4a0 + a0 2 + 2a 2 0) 



T 2 (-2 + 5a - 4a 2 + a 3 ) 
T(-2a + 0-a0 + 3a 2 -a 3 



T (-2a + - a0 + 3a 2 - a 3 ) 
(-20 + 2a0 - 2a 2 + a 3 ) 



(94) 



We see that mn = mn (a, 0, T, N, q 22 ), rfi\i = mi2 (a, 0, T, N, q 22 ) and m 2 2 = ^22 (a, 0, T, N, q 22 ). The usual 
technique for solving the Liapunov equation (84) is by algebraic manipulation and using the symmetry of the matrix, 
as demonstrated with the solution (92). Numerical solutions may be obtained by repeated propagation until steady- 
state is arrived at. 

The time updated steady-state covariance M, referring to (58), is 



We get 



and 



M = lim M (k + l\k) = lim $M (k\k) $' + Q = $ (M N + Mq) $' + Q . 



N 



a(4-2a-0) 
N 



' 1 


T ' 







1 





a(4-2a-0) 



2a 2 + 20- 3a0 (2a - 0) /T 
(2a - 0) /T 20 2 /T 2 

2a 2 + 20 + a0 (2a + 0) /T 
(2a + 0) /T 20 2 /T 2 



1 
T 1 



(95) 
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922 



(-4a/3 + a/3 2 + 2a 2 /3) 



' 1 


T ' 







1 





<?22 



(-4a/3 + a/3 2 + 2a 2 /3) 



T 2 (-2 + 5a - 4a 2 + a 3 ) T (-2a + /3 - a/3 + 3a 2 - a 3 ) 
T(-2a + /3-a/3 + 3a 2 -a 3 ) (-2/3 + 2a/3 - 2a 2 + a 3 ) 

T 2 (-2 + a) T (-2a - /3 + a/3 + a 2 ) 

T (-2a - /3 + a/3 + a 2 ) (-2/3 + 2a/3 - 2a 2 + a 3 ) 



Hence, 



M 



/V 



a (4 - 2a - (3) 



2a 2 + 2/3 + a/3 /3 (2a + /3) /T 
(3 (2a + /3) jT 2/3 2 /T 2 



922 



(-4a/3 + a/3 2 + 2a 2 /3) 
We need steady-state versions of these: designating 



T 2 (-2 + a) T (-2a - /3 + a/3 + a 2 ) 

T (-2a - /3 + a/3 + a 2 ) (-2/3 + 2a/3 - 2a 2 + a 3 ) 





' 





+ 





922 



we then have, referring to (61) and (63), 
then 

Continuing, 
hence, 

Using (81), and (83) 

D = 



D= lim D (k\k) 

k^oo 



D= lim D(k + l\k) 

k— j-oo 



D = FD 
D = D + C 

D=D-C . 

FD=D-C 
FD ~D=(F~I)D 



-C 



D 



1-a (l-a)T " 




' 1 " 




a 




-/3/T 1 - (3 




1 


)"■(- 




) 



-a (l-a)T 
-/3/T -/3 

-/3 -(l-a)T 
/3/T -a 



n -1 r 



a 
/3/T 







a 




' -1 " 










Hence, 



1-a (l-a)T 
-(3/T 1 - (3 





' -1 " 




' a - 1 " 












D = FD = 

Finally, the steady-state time updated total covariance is obtained by substituting into (64), 



Sll S\2 

S21 S22 



Mn M12 
M 2 i M22 



' 1 


T " 




' -1 " 





1 








= M + $T>AT>V 



A-[-l 0] 



1 
T 1 



1 
T 1 



(96) 
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So, 



5 11 S12 

S21 S22 

922 



Mn M12 

M21 M22 

TV 



+ 



A 




a {A -2a- (3) 

n2 



2a 2 + 2/3 + a/3 (3 {2a + (3) /T 
13 (2a + (3) /T 2(3 2 /T 2 



(-4a/3 + a/3 2 + 2a 2 /3) 



T 2 (-2 + a) T (-2a - (3 + a/3 + a 2 ) 

T (-2a - /3 + a/3 + a 2 ) (-2/3 + 2a/3 - 2a 2 + a 3 ) 



In particular, 



" 







' A 


" 





922 


+ 










Sn — 




+ A 





N (2a 2 + 2(3 + a(3) q 22 T 2 (-2 + a) 

a (4 - 2a - /?) (-4a/3 + a/3 2 + 2a 2 /3) 

S21 = M 2 i 

/V/3(2a + /3) q 22 T (-2a - (3 + a[3 + a 2 ) 



a (4 - 2a - (3) T (-4a/3 + a/3 2 + 2a 2 /3) 
Wc turn our attention to (68). In steady-state, 



HS(k + l\k)H' + N = [ 1 



£11 S12 

S2I S 2 2 



N + A 



Sn+N + A 



Also, 



and 



and 



and 



H$D(k\k)A 



[1 0] 



" 1 


T ' 




' -1 " 





1 








• A - 1 • 1 = -A , 



du 



W^- ) AD (k\k)'&H' = -A , 

OA , 



S(k+l\k)H' = 



du 



$D(k\k)A [W 

\ OA 



Sn S12 
S12 S22 



1 T 
1 







Sn 
S12 



A-l-l = 



-A 



Substituting (99), (100), (101), (102) and (103) into (68) gives 

Sn+JV + A-A-A 



a 
(3/T 



Sn 
S12 



+ 



-A 



This, written as two scalar equations, 



a5n+JV + A- A- A =a5n + JV-A =5ii-A 



(97) 



(98) 



(99) 
(100) 
(101) 

(102) 
(103) 
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t ( Su + N + A - A - A) = ^ I 5n + N + -A ) = S 12 ■ 



Substituting (97) and (98) 



/ iV (2a 2 + 2/3 + a/3) q 22 T 2 (-2 + a) \ \ 

[ a(4-2a-0) + (-4a/3 + a/3 2 + 2a 2 (3) J ) 

= N (2a 2 + 2(3 + a/3) q 22 T 2 (-2 + a) A _ A 
a (4 - 2a - (3) (-4a/3 + a/3 2 + 2a 2 /3) 

' /V(2a 2 + 2/3 + a/3) g 22 T 2 (-2 + a) \ ' 

a(4-2a-/3) (-4a/3 + a/3 2 + 2a 2 /3) J , 

_( N (3 (2a + f3) | g 22 T (-2a - /3 + a/3 + a 2 ) \ y 



^ a (4 - 2a - (3) T (-4a/3 + a/3 2 + 2a 2 /3) 
/ N (2a 2 + 2/3 + a/3) g22 r 2 (-2 + a) \ _ N (2a 2 + 2/3 + a/3) 922 T 2 (-2 + a) 



a(4-2a-/3) (-4a/3 + a/3 2 + 2a 2 /3) y a (4 - 2a - /3) (-4a/3 + a/3 2 + 2a 2 /3) 

/ /V (2a 2 + 2(3 + a0) q 22 T 2 (-2 + a) \ = iV/3 (2a + /3) g22 T 2 (-2a - /3 + a/3 + a 2 ) 

^ a(4-2a-/3) (-4a/3 + a/3 2 + 2a 2 /3) J a (4 - 2a - (3) (-4a/3 + a/3 2 + 2a 2 /3) 

We define p = q 22 T 2 /N. With q 22 in (to/ sec) 2 , T in seconds and N in m 2 , p is unitless. Substituting this in the 
previous two equations, we obtain 

^ I (2a 2 + 2(3 + a(3) | p(-2 + a) \_ (2a 2 + 2/3 + a/3) p(-2 + a) 



a(4-2a-/3) (-4a/3 + a/3 2 + 2a 2 /3) J a (4 - 2a - (3) (-Aa(3 + a/3 2 + 2a 2 (3) 

( (2a 2 + 2(3 + a/3) p(-2 + a) \ (3 (2a + (3) p (-2a - (3 + a(3 + a 2 ) 

y a(4-2a-/3) + (-4a/3 + a/3 2 + 2a 2 /3) + y a (4 - 2a - (3) + (-4a/3 + a/3 2 + 2a 2 /3) ' 

We divide the previous two equations and cancel an a, 

a _ (2a 2 + 2/3 + a/3) (-4/3 + /3 2 + 2a/3) + p (-2 + a) (4 - 2a - /?) 

/3 _ /3 (2a + (3) (-4/3 + /3 2 + 2a/3) + p (-2a - /3 + a/3 + a 2 ) (4 - 2a - /3) ' 

Cross multiplying gives: 



2/3 4 + (4a - 8) (3 3 + p ((a 2 - 2a + 2) (3 2 + (3a 3 - 10a 2 + 12a - 8) (3 + (2a 4 - 8a 3 + 8a 2 )) = . (104) 

Equation (104) gives our relationship between a and (3. The noise ratio p is a parameter in the equation which is 
known, or at least known to be within a range. Equation (104) may be factored 

((3 + (2a - 4)) (2/3 3 + p ((a 2 - 2a + 2) (3 + a 2 (a - 2))) = . (105) 

Hence, (3 = (4 — 2a) is a solution that is independent of p. This solution is not permitted however since it 
violates condition 3. of (91). There is a second real solution for (3 given a (which depends on p.) The remaining 
two solutions for (3 given a may be a complex conjugate pair. The table below gives some representative solutions 
to (104). The two real solutions are presented. The second one listed we don't use because of the condition 3. The 
Newton- Raphson method may be used to compute all of the solutions to (104). 
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p 



a 



2 
4 
6 
6 
8 
8 
10 
10 
10 



0.2 
0.2 
0.2 
0.4 
0.2 
0.4 
0.2 
0.4 
0.5 



0.04385, 3.6 
0.04386, 3.6 
0.04389, 3.6 
0.1866, 3.2 
0.04389, 3.6 
0.1870, 3.2 
0.04389, 3.6 
0.1873, 3.2 
0.2959, 3.0 



These a and f3 give that the eigenvalues of F, in (82), have norm less than 1. Hence, by Theorem 2.1, page 64 
of [1], Mn and Mq, the solutions to (86) and (87) respectively, exist are unique and are positive definite. 

4 Summary and Conclusions 

In this paper we considered some topics in radar sensor bias. We presented an algorithm that estimates the absolute 
bias of two sensors when the relative bias between the sensors is given. The algorithm uses the relative bias, which 
is given in rectangular coordinates, as a constraint. The absolute biases, in spherical coordinates, for the sensors are 
obtained by the solution to an optimization problem that exploits the spherical-to-rectangular coordinate conversion. 
We presented a reduced-state filter that is designed for performance with sensor bias. The filter is reduced-state 
since it does not contain additional bias states. The filter design is influenced by the filter in [4]. It may be viewed 
as a dual design (in the control theory sense) to the filter in [4] . 

A flow diagram for processing radar data with bias may contain these stages: 

1. Estimate state with the a — (3 filter optimized for measurement bias, as presented in Section 3. 

2. For a multi-sensor problem, estimate the relative sensor bias using an optimized algorithm such as in [2]. 

3. Continue by estimating the absolute bias for each sensor using the algorithm presented in Section 2. 
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In this appendix we present the transformation from the ENU{1) to the ENU(2) coordinate systems. But first, 
consider the transformation from ECI to ENU. Consider an ENU coordinate axis located at longitude-latitude 



1997. 
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Appendix: Transformation from ENU(1) to ENU{2) 
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O — L and define the rotation matrix 
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EC I 



\J\ — e 2 sin 2 L 



Tenu2ECiPenu 



with Tenu2ECI = T E ci2ENU- For position, we need to include a translation, so that for a given position vector in 
ENU coordinates 

cos L cos 
cos L sin 
(1 - e 2 ) sini 

where r ee is the earth's equatorial radius and e is the earth's eccentricity. For velocity, use the rotation alone. 

Let the ENU(1), ENU(2) coordinate system be located at longitude-latitude Oi — L\ and f^2 — L2 respectively. 
Next we consider our transformation going from Oi — L\ to O2 — L2- The rotation part of this transformation can be 
represented by the matrix below (going from ENU (1) to ENU {2)). There are 3 steps. First, the ENU (1) coordinates 
are rotated down to the equator. Second, these coordinates are rotated along the equator by the longitude difference. 
Third is the rotation up to the latitude of the ENU (2) system. 



-ENU(l)2ENU(2) 



1 

cos L2 — sin L 2 
sin L 2 cos L 2 



cos (0 2 ~ f^i) 



— sin (O2 


cos (O2 — ^1) 



cos (O2 — ^1) 
— sin_L2 sin (O2 — Oi) 
cos L2 sin (0 2 — 1 ) 



(Note 





1 
sin(fi 2 -fii) 

sinLi sin (0 2 — £li) 
cos L\ cos Z/2 + sinLi sin L2 cos (fl 2 — 0\) cos L2 sin L\ 
cosLi sinL 2 — cosL 2 sinLi cos (0 2 — 0\) sinLi sinL 2 



1 

cosLi sinLi 
— sinii cosLi 



— cos L\ sin (Q2 — Oi ) 
in Li — cos L\ sin L2 cos (O2 — Oi) 
; ■■■ ' Li cosL 2 cos (f2 2 — ill) 



- cos. 



cos (0 2 — fli) — sin (O2 — fii) 

1 

sin(f22 — Oi) cos(^2 — ^1) 



cosfii cosfl 2 + sin ill sinfl 2 — cosfii sinil 2 + cos£! 2 sinOi 

1 

cos lli sin Sl 2 — cos f2 2 sinfii cos Qi cos fl 2 + sin Oi sin 2 



cosJli sinfli 

1 
— sin fii cos fli 



cos SI2 — sin fl2 

1 
sin O2 cos V.j 

so that the rotation T E nu(i)2ENU(2) is a rotation from the first coordinates down to 
coordinates.) We have 

TeNU(2)2ENU(1) = T' E NU(1)2ENU{2) 

The position vector from the ENU(1) to the ENU(2) coordinate axes (in ECI coordinates) is 

cos L2 cos 0,2 



ECI and then up to the second 



and 



PeNU(1)2ENU(2),ECI 

in the other coordinates this vector is 



— ^ee ' 



cos L 2 sin O2 
(l - e 2 ) sinL 2 



cos L\ cos 0\ 
cosLi sin Hi 
(l - e 2 ) sinLi 



/ o 2 ee / 2 

V 1 — e 2 sin L2 \/\ — e 2 sin L\ 

PENU(l)2ENU(2),ENU(i) = T E CI2ENU(i)PENU(l)2ENU(2),ECI 



ition coordinate transformation, including translation can be represented by 



The total position 

PeNU(2) — —PeNU(1)2ENU(2),ENU(2) + T E NU(\)2ENU(2)PeNU(\ 

The total velocity coordinate transformation is given by the rotation alone. 
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